{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 8\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions from the previous chapter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_results(census, un, timeseries, title):\n",
    "    \"\"\"Plot the estimates and the model.\n",
    "    \n",
    "    census: TimeSeries of population estimates\n",
    "    un: TimeSeries of population estimates\n",
    "    timeseries: TimeSeries of simulation results\n",
    "    title: string\n",
    "    \"\"\"\n",
    "    plot(census, ':', label='US Census')\n",
    "    plot(un, '--', label='UN DESA')\n",
    "    plot(timeseries, color='gray', label='model')\n",
    "    \n",
    "    decorate(xlabel='Year', \n",
    "             ylabel='World population (billion)',\n",
    "             title=title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Simulate the system using any update function.\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that computes the population next year\n",
    "    \n",
    "    returns: TimeSeries\n",
    "    \"\"\"\n",
    "    results = TimeSeries()\n",
    "    results[system.t_0] = system.p_0\n",
    "    \n",
    "    for t in linrange(system.t_0, system.t_end):\n",
    "        results[t+1] = update_func(results[t], t, system)\n",
    "        \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_table2(filename = 'data/World_population_estimates.html'):\n",
    "    tables = pd.read_html(filename, header=0, index_col=0, decimal='M')\n",
    "    table2 = tables[2]\n",
    "    table2.columns = ['census', 'prb', 'un', 'maddison', \n",
    "                  'hyde', 'tanton', 'biraben', 'mj', \n",
    "                  'thomlinson', 'durand', 'clark']\n",
    "    return table2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#table2 = read_table2()\n",
    "#table2.to_csv('data/World_population_estimates2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "table2 = pd.read_csv('data/World_population_estimates2.csv')\n",
    "table2.index = table2.Year\n",
    "table2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "un = table2.un / 1e9\n",
    "census = table2.census / 1e9\n",
    "plot(census, ':', label='US Census')\n",
    "plot(un, '--', label='UN DESA')\n",
    "    \n",
    "decorate(xlabel='Year', \n",
    "             ylabel='World population (billion)',\n",
    "             title='Estimated world population')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the quadratic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the update function for the quadratic growth model with parameters `alpha` and `beta`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func_quad(pop, t, system):\n",
    "    \"\"\"Update population based on a quadratic model.\n",
    "    \n",
    "    pop: current population in billions\n",
    "    t: what year it is\n",
    "    system: system object with model parameters\n",
    "    \"\"\"\n",
    "    net_growth = system.alpha * pop + system.beta * pop**2\n",
    "    return pop + net_growth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the starting time and population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = get_first_label(census)\n",
    "t_end = get_last_label(census)\n",
    "p_0 = get_first_value(census)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the system object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "system = System(t_0=t_0, \n",
    "                t_end=t_end,\n",
    "                p_0=p_0,\n",
    "                alpha=0.025,\n",
    "                beta=-0.0018)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the model and plot results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, update_func_quad)\n",
    "plot_results(census, un, results, 'Quadratic model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate projections, all we have to do is change `t_end`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "system.t_end = 2250\n",
    "results = run_simulation(system, update_func_quad)\n",
    "plot_results(census, un, results, 'World population projection')\n",
    "savefig('figs/chap08-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The population in the model converges on the equilibrium population, `-alpha/beta`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[system.t_end]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "-system.alpha / system.beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  What happens if we start with an initial population above the carrying capacity, like 20 billion?  Run the model with initial populations between 1 and 20 billion, and plot the results on the same axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the projection from our model with projections produced by people who know what they are doing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_table3(filename = 'data/World_population_estimates.html'):\n",
    "    tables = pd.read_html(filename, header=0, index_col=0, decimal='M')\n",
    "    table3 = tables[3]\n",
    "    table3.columns = ['census', 'prb', 'un']\n",
    "    return table3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#table3 = read_table3()\n",
    "#table3.to_csv('data/World_population_estimates3.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "table3 = pd.read_csv('data/World_population_estimates3.csv')\n",
    "table3.index = table3.Year\n",
    "table3.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`NaN` is a special value that represents missing data, in this case because some agencies did not publish projections for some years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function plots projections from the UN DESA and U.S. Census.  It uses `dropna` to remove the `NaN` values from each series before plotting it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_projections(table):\n",
    "    \"\"\"Plot world population projections.\n",
    "    \n",
    "    table: DataFrame with columns 'un' and 'census'\n",
    "    \"\"\"\n",
    "    census_proj = table.census / 1e9\n",
    "    un_proj = table.un / 1e9\n",
    "    \n",
    "    plot(census_proj.dropna(), ':', color='C0', label='US Census')\n",
    "    plot(un_proj.dropna(), '--', color='C1', label='UN DESA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the model until 2100, which is as far as the other projections go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = System(t_0=t_0, \n",
    "                t_end=2100,\n",
    "                p_0=p_0,\n",
    "                alpha=0.025,\n",
    "                beta=-0.0018)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, update_func_quad)\n",
    "\n",
    "plt.axvspan(1950, 2016, color='C0', alpha=0.05)\n",
    "plot_results(census, un, results, 'World population projections')\n",
    "plot_projections(table3)\n",
    "savefig('figs/chap08-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "People who know what they are doing expect the growth rate to decline more sharply than our model projects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:** The net growth rate of world population has been declining for several decades.  That observation suggests one more way to generate projections, by extrapolating observed changes in growth rate.\n",
    "\n",
    "The `modsim` library provides a function, `compute_rel_diff`, that computes relative differences of the elements in a sequence.\n",
    "\n",
    "Here's how we can use it to compute the relative differences in the `census` and `un` estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_census = compute_rel_diff(census)\n",
    "plot(alpha_census, label='US Census')\n",
    "\n",
    "alpha_un = compute_rel_diff(un)\n",
    "plot(alpha_un, label='UN DESA')\n",
    "\n",
    "decorate(xlabel='Year', label='Net growth rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other than a bump around 1990, net growth rate has been declining roughly linearly since 1965.  As an exercise, you can use this data to make a projection of world population until 2100.\n",
    "\n",
    "1. Define a function, `alpha_func`, that takes `t` as a parameter and returns an estimate of the net growth rate at time `t`, based on a linear function `alpha = intercept + slope * t`.  Choose values of `slope` and `intercept` to fit the observed net growth rates since 1965.\n",
    "\n",
    "2. Call your function with a range of `ts` from 1960 to 2020 and plot the results.\n",
    "\n",
    "3. Create a `System` object that includes `alpha_func` as a system variable.\n",
    "\n",
    "4. Define an update function that uses `alpha_func` to compute the net growth rate at the given time `t`.\n",
    "\n",
    "5. Test your update function with `t_0 = 1960` and `p_0 = census[t_0]`.\n",
    "\n",
    "6. Run a simulation from 1960 to 2100 with your update function, and plot the results.\n",
    "\n",
    "7. Compare your projections with those from the US Census and UN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Related viewing:** You might be interested in this [video by Hans Rosling about the demographic changes we expect in this century](https://www.youtube.com/watch?v=ezVk1ahRF78)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
